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TECHNICAL PAPER 


AN EFFICIENT ALGORITHM FOR GENERATING RANDOM NUMBER PAIRS 
DRAWN FROM A BIVARIATE NORMAL DISTRIBUTION 

INTRODUCTION 


An efficient algorithm was devised to generate random numbers drawn from a bivariate normal dis- 
tribution with any desired values of the two means, two standard deviations, and correlation coefficient. 
A sequence of pairs of uniformly distributed random numbers are transformed to obtain a sequence of 
random number pairs which are bivariate normally distributed. The only accuracy limitations are in the 
quality of the uniform random number generation routine, and in the ability of the computer to perform 
exact calculations. 

The present technique is a special case of a general method described by Scheuer and Stoller [ 1 ] 
who devised a method for generating sequences of random n-vectors drawn from a multivariate normal 
distribution with any desired covariance matrix. For an n-vector a sequence of n independent normal variates 
was generated. A linear transformation was applied to each of the n-vectors of the sequence, and the result 
was a sequence with the desired properties. Since the transformation was linear the transformed sequence 
was also linear, but the vector components are no longer uncorrelated. By judicious selection of the trans- 
formation matrix, the desired covariance matrix was obtained. For the present case, the linear transforma- 
tion involves multiplication by scale factors and rotation. 


MATHEMATICAL DEVELOPMENT 

The bivariate normal density function with zero means is given by: 

x^ 2rxy y^ 
a x 2 CT x a y a 2 

where 

o 2 = E(x 2 ) 
a y 2 = E(y 2 ) 
r = E(xy)/a x a y 


^x,y 


2tt ct x a yA /l-r 2 


exp 


' 2\/l-r2 


( 1 ) 


and E( ) means the expected value of the variable in parenthesis. The means can be set to zero without loss 
of generality because transforming to variables with nonzero means involves a simple addition. 



The basis for the algorithm is the generation of a sequence of independent normally distributed 
random pairs (x',y') (being independent the correlation coefficient r is equal to zero) which are transformed 
to obtain a new sequence (x,y) which is normally distributed with the desired standard deviations and 
correlation coefficients. For nonzero means, the desired means are then added to each pair (x,y). Because 
rotation, multiplication by scale factors, and addition constitute a linear transformation, the resulting pairs 
will be normally distributed. Any desired method of generating (x',y') can be used, but for convenience, a 
method described by Box and Muller [2] was utilized. Other methods involving the central limit theorem 
and rational approximation would be acceptable. Some of these techniques are described by Howell and 
Rheinfurth [3], Rheinfurth [4] has also pointed out that the Box and Muller method may have some 
nonzero values of serial correlation possibly resulting from function evaluation errors or other computer 
arithmetic inaccuracy. In any event any normally accepted method for generating random normal deviates 
should work for the current algorithm. 

The technique of Box and Muller generates a pair of independent normal deviates by means of the 
following transformation: 


x" = (- 21 nuj )'*/2 s j n 2?r U 2 
y" = (-21 mi})" cos 2 n U2 


( 2 ) 


u j and U 2 are a pair of random variates uniformly distributed between 0 and 1. The inverse transforma- 
tions are given by: 


uj = exp[-(x" 2 + y" 2 )/2] 


U 2 = arc tan (y" / x”)/( 2 tt') 


(3) 


The Jacobian of these inverse transformations is: 

-exp[-(x" 2 + y" 2 )/2] 


J[(ui,u 2 )/(x",y )] = 


2n 


(4) 


The joint distribution of x" and y" is given by: 


f x",y" 


(x",y") = f Ul u 2 Jt( u 1> U 2 >/( x ",y")] 


(5) 


or 


f x",y"(x",y") = l/(>/5o exp(-x" 2 /2) \/(^) exp(-y" 2 /2) 


( 6 ) 
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From the above equation it is obvious that x" and y" are independent normally distributed variables with 
zero means and unit standard deviations. 


The transformations in theory are exact. In practice computer errors of arithmetic, function evalua- 
tion, etc. create some inaccuracies which are believed to be small. 

Equation (2) generates a sequence of normally distributed random number pairs. Transformation to 
the desired sequence requires three simple operations. The first is the multiplication given in the following 
equation : 


x' = a x ' x ” 

y' = <y y" 


(7) 


The standard deviations of x' and y' are a x > and ay, respectively. The actual values of o x * and ay are yet 
to be determined. 

The next transformation is a rotation. Figure 1 illustrates the problem, x and y are given by the 
rotational transformation: 


x = x' cos 0 - y' sin 0 

(8a) 

y = x' sin 8 + y' cos 8 

(8b) 



Figure 1. Rotation of axes transformation. 
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The inverse transformation is: 


x' = x cos 9 + y sin 9 


(9a) 


y' = -x sin 9 + y cos 9 


(9b) 


Multiplying (9a) by (9b) and taking expected values gives: 

a x ' <Ty' r' = (a y ^ - a x 2 ) cos 0 sin 0 + a x Oy r (cos 0 - sin 0) . (10) 

r' is the rotated correlation coefficient and is therefore zero. With the left hand side of (10) set to zero the 
equation can be solved for 9 in terms of the known quantities a x , a y , and r. The result is: 

0 =(l/ 2 )arctan [ 2 ra x a y /(CT x 2 -ff y 2 )] . (11) 

Now that 9 is known, expressions for the unknown quantities a x > and a y ' in equation (7) are required. 
Multiplying equations (9a) and (9b) together and taking expected values yields the following results: 

°x' 2 = °x‘ cos ^ 9 + 2r o x Oy cos 9 sin 9 + o y 2 sin 2 9 (12a) 

0 y ' 2 = a x 2 sin 2 9 - 2r a x o y cos 9 sin 9 + a y 2 cos 2 9 . (12b) 

Notice that o x > and Oy> are now expressed in terms of known quantities and can be evaluated. Adding (12a) 
to ( 1 2b) gives the following relation : 

o x ' 2 + Oy' 2 = c x 2 + Oy 2 . (13) 


If double angle formulae are substituted into (12a) and (12b) and the equations differentiated with respect 
to 9 and the resulting derivatives set to zero, for both (12a) and (12b), equation (11) results. This means 
that for the standard deviations a x > and a y ' of the rotated variates, extrema occur when x' and y' are inde- 
pendent. If x' corresponds to the probability major axis and y' to the minor axis then the value o x < is the 
maximum standard deviation for any rotation angle and Oy> is a minimum. 

From equation (11), 9 is known, and from two of the three equations (12a), (12b), and (13), a x > 
and Oy> are known. Each of the three parameters a x q a y », and 9 are known in terms of the desired parame- 
ters a x , Oy, and r. With this information, (x,y) can be calculated. 

All required equations are now at hand, and the desired sequence can be generated. The algorithm 
for generating the desired variate pairs is summarized below. 
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1 ) 

program. 

2 ) 

3) 

4) 

5) 

6 ) 

7) 

8 ) 
9) 


Input desired value of the two means, two standard deviations, and correlation coefficient to the 

Use these values to calculate the rotation angle, using (11). 

Calculate ct x ' and Oy> using equations (12a) and (12b). 

Generate two random numbers (uniformly distributed) Uj, and U2* 

Calculate x", and y" from equation (2). 

Calculate x', and y' from equation (7). 

Calculate x and y from equation (8a) and (8b). 

Add the desired mean values to x and y„ 

If another random number pair (x,y) is needed return to step 4; otherwise stop. 

TEST OF THE ALGORITHM 


The previously described algorithm was tested with the computer program described in the 
Appendix. The code was written for a Hewlett-Packard 1000 F Series computer. Results are summarized 
in Tables 1 and 2. Table 1 compares the desired values of o x , a y, and r with those calculated for the 

sequence. Each sequence consisted of 1000 number pairs. The error is quite small. The purpose of Table 2 
is to show that as the size of the sample increases, the deviations from desired values decrease. 

TABLE 1. COMPARISON OF STATISTICS OF BIVARIATE PAIRS WITH THE 
DESIRED PARAMETER VALUES (N = 1000 FOR EACH CASE) 


°x 

°xm 

°y 

ym 

i 

T 

m 

Ar 

1.25 

1.2379 

1.0 

1.0002 

0.25 

0.2538 

0.0038 

1.25 

1.2376 

1.0 

0.9996 

0.50 

0.5000 

0.0000 

1.25 

1.2382 

1.0 

0.9979 

0.75 

0.7487 

-.0013 

1.50 

1.4871 

1.0 

0.9998 

0.25 

0.2562 

0.0062 

1.50 

1.4863 

1.0 

1.0001 

0.50 

0.5018 

0.0018 

1.50 

1.4867 

1.0 

0.9984 

0.75 

0.7494 

-.0006 

2.00 

1.9845 

1.0 

0.9993 

0.25 

0.2576 

0.0076 

2.00 

1.9836 

1.0 

1.0002 

0.50 

0.5035 

0.0035 

2.00 

1.9837 

1.0 

0.9990 

0.75 

0.7501 

0.0001 

3.00 

2.9782 

1.0 

0.9990 

0.25 

0.2582 

0.0082 

3.00 

2.9775 

1.0 

1.0002 

0.50 

0.5045 

0.0045 

3.00 

2.9774 

1.0 

0.9994 

0.75 

0.7507 

0.0007 

4.00 

3.9714 

1.0 

0.9989 

0.25 

0.2584 

0.0084 

4.00 

3.9710 

1.0 

1.0002 

0.50 

0.5047 

0.0047 

4.00 

3.9708 

1.0 

0.9996 

0.75 

0.7509 

0.0009 

5.00 

4.9646 

1.0 

0.9988 

0.25 

0.2584 

0.0084 

5.00 

4.9642 

1.0 

1.0002 

0.50 

0.5049 

0.0049 

5.00 

4.9641 

1.0 

0.9996 

0.75 

0.7510 

0.0010 

10.0 

9.9300 

1.0 

0.9988 

0.25 

0.2585 

0.0085 

10.0 

9.9298 

1.0 

1.0001 

0.50 

0.5050 

0.0050 

10.0 

9.9298 

1.0 

0.9997 

0.75 

0.7512 

0.0012 
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TABLE 2. CONVERGENCE OF STATISTICS TO DESIRED VALUES FOR 
INCREASING NUMBERS OF VARIATES 


N 

Aa x 

A a 

y 

10 

0.7962 

0.7700 

100 

-0.038 

0.0794 

1000 

-.0164 

0.0002 

10000 

-.0096 

0.0044 


Ar 


0.1787 

0.1267 

0.0035 

0.0023 


Figures 2 and 3 show interesting trends in the error of the correlation coefficient. In Figure 2 corre- 
lation coefficient error Ar, which is the difference in the calculated correlation and the desired correlation, 
is plotted versus the desired correlation coefficient. The error decreases from a maximum at r = 0 to near 
zero at r = ±1. This is for the case a x = 2, Oy = 1. Figure 3 plots the correlation coefficient error versus 

values of the ratio o x /o y. Each curve corresponds to a different value of the desired correlation coefficient, 

r. The regularity of the curves depicted in Figures 2 and 3 was somewhat surprising at first glance, but a 
relatively simple model was developed which explained the qualitative aspects of the plots. 

Returning to the transformations (9a) and (9b), the value of the rotated correlation coefficient, r„ 
is obtained by multiplying the equations together, taking expected values, and manipulating. 


sin 20 (o ^ _ a x ^)l2 + a x a y cos 26 

r' = . (14) 

<J V ' O.J 


Values of o x > and Oyt are obtained from equations (12a) and (12b). The variation of r' with o x , <jy, 6, and 

r is shown in Figures 4 and 5. The value of the unrotated correlation coefficient, r, for each of the curves is 
r' evaluated at 6 = 0. For each curve four zeroes of r' are seen. The value of 6 corresponding to the zero 
closest to 6 = 0 is the value of 6 used in transformations (9a) and (9b). For a positive value of the correla- 
tion coefficient, r, the closest zero corresponds to a positive value of 6. For a negative value of r, the 
corresponding value of 6 is negative. The four zeroes of r' are expected since the x' axis will correspond to 
the positive and negative major and minor axes of the probability ellipse for four different values of 6. 

Looking at the curves of Figures 4 and 5 and paying particular attention to the zero-crossing behavior 
of r\ one might expect that the error in r may be related to the value of the derivative of r' at the origin. 
The x'-y' axes are rotated by the angle to obtain the desired values of o x , a y, and r. Some errors are involved 

in the transformation and an error in 6 results in an error in r; i.e., r = dr'/d0 Ad . In this model Ar varies in 
magnitude with dr'/d0 evaluated at the origin. 

Differentiating equation (14) with respect to 6 and evaluating at 6 = 0 gives the following result: 


dr 

d0 


0=0 


= (°y/ a x - a x/°y) 0 ~r 2 ) 


(15) 


If the ratio a x /a y is held constant, (15) can explain qualitatively the behavior depicted in Figure 2. The 
maximum error occurs at r = 0, and minimum at r = ±1 as shown. 
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Figure 3 is not totally explained by the model, however. The sharp decrease of Ar near o x /a y = 0 
is in good agreement with the model, but as <7 x /oy 00 the curves seem to level off approaching a constant 

value. Equation (15) predicts that the error increases linearly with a x /oy. The reason for the discrepancy is 
not known, but the actual behavior is better than predicted. 



Figure 2. Variation of the absolute error of the correlation coefficient (a x = 2, a y = 1). 


r = 0.25 


Ar 



© © .0051 


© Q .0086 


A .0012 


1000. 


Figure 3. Variation of correlation error with ratio of standard deviations. 
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Figure 4. Variation of rotated correlation coefficient with standard deviation ratio 

(o x /a y = 1,2,5, 100; r= 0.5). 
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Figure 5. Variation of unrotated correlation coefficient with rotation angle and unrotated 
correlation coefficient (a x /ay = 2, r = 0.0, 0.5, 0.9, 1.0). 


Qualitatively the increase in Ar with decreasing r as shown in Figure 3 is in agreement with (15). 
The ratios of corresponding values Ar (r = 0.25); Ar (r = 0.5); Ar (r = 0.75) are not in quantitative agree- 
ment, however. 


SUMMARY 


An algorithm for drawing random number pairs from a bivariate normal distribution with desired 
means, standard deviations, and correlation coefficient was described. The algorithm is efficient in the sense 
that a sequence of uniformly distributed random number pairs is converted into a sequence of the same 
number of bivariate normally distributed random number pairs. In other words, for each pair of random 
numbers the algorithm generates a pair of bivariate normal random numbers. 

A FORTRAN program was written to test the algorithm. The algorithm worked well, with some 
small errors in the parameters. The regularity of the errors seems to discount the possibility that the errors 
are entirely statistical in nature. 

To better understand the errors, a simple model was devised which qualitatively explained the 
error in the correlation coefficient. The model was in qualitative agreement with most aspects of the error, 
and when the model was not in agreement, the actual behavior was better than predicted. 
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APPENDIX 


FORTRAN PROGRAM FOR GENERATING BIVARIATE 
NORMAL RANDOM NUMBERS 


The listing is of a FORTRAN program which generates bivariate normal random numbers. The 
routine RMPAR is an HP RTE routine which is described in Reference 5. RMPAR passes as many as five 
integer parameters to the program through the array IPAR. The calling sequence is 

RU,BIVAR,IPRT1 5 IPRT2,NRV 

where IPRT1 is the list device where the statistical summary of the random number sequence is sent IPRT2 
is the list device where the random variates are sent, and NRV is the number of variate pairs generated. The 
program is written for interactive input from a CRT terminal which is logical unit number 1. 

System routine URAN generates uniformly distributed random numbers, and is described in Refer- 
ence 6 0 : 


FTN4X, L 

PROGRAM B I VAR 


* 

A 

* 

* 


* 


:+i 


* 

* 

* 

* 

* 


* A A A +• * A * A A A +. A A A A A A A * A A'- A A A: A A :+: A A A- i ■ - t- t= A- -1 A A A- A t- A A A Y A A -A '+ -A A A -A ’A ■ A A- -A -A 

■A -A A * 

** * THIS PROGRAM GENERATES BIVARIATE NORMAL RANDOM NUMBERS. * 

A -A A * 

a A a PROGRAMMER: NAP REN CAMPBELL * 

:A A *' 

A -A A c A L L I H G S E Q IJ t£ N C E : R U . h I V A k , I P k T 1 . 1 P h* T £ . H R v * 

* + * V A R I A B L E D F. F I H I T I 0 H ; + 


** * SIGN = DESIRED X STD DEVIATION 

** * SI GY = DESIRED V STD DEVIATION 

** a R = DESIRED CORRELATION COEFFICIENT 

** * THETA = ANGLE TO ROTATE AXES 

a a a X = X VALUE OF BIVARIATE NORMAL RANDOM HUMBER 

A a a V = Y VALUE OF BIVARIATE NORMAL RANDOM NUMBER 

aa a XPR = ROTATED X VALUE < I NDEP OF YPR > 

a a a YPR = ROTATED V VALUE < I NDEP OF XPR > 

aa a IPRT1 = STATISTICS SUMMARY PRINTER NUMBER 

AA A IPRT2 = NORMAL VARIATES PRINTER NUMBER 

AA a NRV = NUMBER OF RANDOM VARIATES TO BE GENERATED 

-+: A A 


A 


A 

A 


A 

A 

A 


A 


A 


■A A A A A A A A % A A A A A A A A A A A A A A A A A A * A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A 


A A 


A A 


DIMENSION I PARC 5) 
CALL RMPAR** IPAR) 

A A A A A A A A A A AAA A A A A A 

IPRT1 = I PAR< 1 > 
IPRT2 = IPAR*: 2> 

NRV = I PARC 3) 

A A A A A A A A A A A A A A A A A A 
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TWQPI = 2 . *3. 1 41 592654 
SUMX = 0. 

SUMXX = 0, 

SUMY = 0. 

SUMYY =0. 

SUMXY = 0. 

WRITE-;' 1 , 9999 > 

9999 FORMAT < " ENTER DES I RED S I GX , S I GY , R$ " > 

RE ADC 1 , * ) SIGX, SI GY , R 
WRITE< I PRT 1 , 9998 > 

99 98 F ORMAT C //// > 

WRITE< IPRT1 , 9990 > 

99 9 0 F ORMAT -; / " * **** + +-<< + + + + ♦+ + * + * :f. + + * * + •+ ■+ + + ■+: + + + '+ :*-■ * + 


WRITE-; IPRT1 ,9991 > NRV 

9991 FORMAT < " NUMBER OF BIVARIATE PAIRS GENERATED IS “15) 

THETA = 0 , 5+AT ANC 2 . *R*S I GX+S I GY/C S I G X * 2 - S I G Y + *2 > ) 

A = SIN-;' THETA) 

B = COS< THETA) 

DTHTA = 360 . *THETA/TWOPI 
WR I TE< I PRT 1 , 9997 ) S I GX , S I GY , R , DTHTA 
9997 FORMAT< " SIGX = "F8 . 3“ SIGY = "F8.3” R =- "F7.4" THETA = " 
* F 9 , 4 " DEGREES" > 

WRITE-; I PRT 1 , 9989) 

9989 FORMATC 7/ > 

SXPR = <B*SIGX>**2 + 2 , * R * S I G X * S I G V * A + B + CSIGY*A)**2 
SXPR = SQRT( SXPR ) 

SYPR = < A*S I GX >**2 ~ 2 . * R * S I GX*S I GY+A+B + C S I GY*B > +• * 2 
SYPR = SORT-; SYPR) 

WR I TE< I PRT 1 , 9996 > SXPR , SYPR 
9996 FORMATC " SXPR = "F9.4" SYPR = "F9,4> 

WRITE< IPRT1 , 9989 > 

WRITE< IPRT2 , 9995 ) 

9995 FORMATC “ I X Y SUMXX 

+ SUMYY SUMXY" ) 

DO 100 1=1, NRV 
U1 = UR AN-; 1 ) 

U2 = URAN< 1 >*TWGPI 

TEMP = S 0 R T < - 2 , * A L 0 G C U t ) ) 

XPR = TEHP-+SIHC U2 )*SXPR 
Y'PR = TEMP*COSC U2 >*SYPR 
X = XPR+B - YPR*A 
Y = XPR*A + Y P R * B 
SUMX = SUMX + X 
SUMY = SUMY + Y 
SUMXX = SUMXX + X+X 
SUMYY = SUMYY + V*Y 
SUMXY = SUMXY + X*Y 

WR I T E < I P R T 2 , 9 9 9 4 > I , X , Y , S U M X X , S U M Y Y , S U M X Y 
9994 FORMATC 5X , 1 5 , 5C 5X , F 1 0 , 5 ) ) 

100 CONTINUE 
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XBAR = SUMX/NRV 
YBAR = SUMY/NRV 
STDX = SQRT* SUMXX/* HRV - 1 . > > 

STDY = SQRT* SUMYY/* HRV - 1 , > > 

RM = SUMXY/<NRV*STDX*STDY> 

WR I TE< IPRTt , 9989 > 

WRITE* IFRT1 , 9993) YBAR, YBAR 
9993 FORMAT* " XBAR = "F8,4 n YBAR = ”F8.4> 

WR I TEC I PRT I j 9992 > S I GX , S I GY , R , STDX , STDY , RM 
9992 FORMAT*” SI GX = ”F9 , 4 11 SIGY= W F9 . 4 ” R= i! F7,4 

* F9 , 4 ” RM="F7.4> 

WRITE* I PRT 1 , 9990 > 

WRITE* IPRT2,9998> 

STOP 

END 

END* 


STDX- !I F9 , 4" STDY= u 
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